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Abstract 

I present a numerical package (CosmoTransitions) for analyzing finite-temperature 
cosmological phase transitions driven by single or multiple scalar fields. The 
package analyzes the different vacua of a theory to determine their critical 
temperatures (where the vacuum energy levels are degenerate), their super- 
cooling temperatures, and the bubble wall profiles which separate the phases 
and describe their tunneling dynamics. 1 introduce a new method of path 
deformation to find the profiles of both thin- and thick-walled bubbles. Cos- 
moTransitions is freely available for public use. 

Keywords: Finite-temperature Field Theory, Cosmology, Phase 
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OS with Python installed. 

RAM: Approximately 50 MB, mostly for loading plotting packages. 
Keywords: Finite-temperature Field Theory, Cosmology, Phase Transitions 
Classification: 1.9, 11.1 

External routines/libraries: SciPy, NumPy, matplotLib 

Nature of problem: I describe a program to analyze early- Universe finite-temperature 
phase transitions with multiple scalar fields. The goal is to analyze the phase 
structure of an input theory, determine the amount of supercooling at each phase 
transition, and find the bubble-wall profiles of the nucleated bubbles that drive 
the transitions. 

Solution method: To find the bubble-wall profile, the program assumes that tun- 
neling happens along a fixed path in field space. This reduces the equations of 
motion to one dimension, which can then be solved using the overshoot/undershoot 
method. The path iteratively deforms in the direction opposite the forces perpen- 
dicular to the path until the perpendicular forces vanish (or become very small). 
To find the phase structure, the program finds and integrates the change in a 
phase's minimum with respect to temperature. 

Running time: Approximately 1 minute for full analysis of the two-dimensional 
test model on a 2.5 GHz CPU. 



1. Introduction 

Phase transitions driven by scalar fields likely played an important role 
in the very early evolution of the Universe. In most inflationary models, 
the dynamics are driven by the evolution of a scalar infiaton field, while 
at later times electroweak symmetry breaking is thought to be driven by a 
transition in the Higgs field vacuum expectation value. Electroweak scale 
physics is currently being probed by the LHC, so the phenomenology of the 
electroweak phase transition is of particular interest. A strongly first-order 
electroweak phase transition would have been a source of entropy production 
in the early Universe (thereby changing the evolution of its scale with re- 
spect to temperature) and produced a stochastic background of gravitational 
radiationJT], perhaps observable by future space-based gravitational radiation 
observatories [2]. In addition, a strongly first-order electroweak phase tran- 
sition may have satisfied the Sakharov conditions [3] and been responsible 
for the current baryon asymmetry of the universe (for recent studies see e.g. 
Refs. HElElIIlElElIinilllliniliaiiailSl [IS]), or may have affected the 
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relic density of (for example) dark matter particles [T71 HSl I 

In the standard model, the electroweak phase transition is not strongly 
first-order unless the Higgs mass is below ~ 70 GeV[20l EU |22], which is 
excluded by the current LEP bound of 114.4 GeV[23]. However, electroweak 
baryogensis can be saved in extensions to the standard model, many of which 
include extra dynamic scalar fields (such as two-Higgs-doublet models [211 
l25l [26l ETJ EHl EHl [30]). The amount of produced baryon asymmetry depends 
crucially upon the dynamics of the phase transition, and particularly upon 
the bubble-wall profiles that separate the high- and low-temperature phases. 
These profiles are fairly easy to calculate if there is only one scalar field, but 
multiple fields greatly increase the computational complexity. 

In this paper, I present an easy-to-use numerical package (CosmoTran- 
sitions) to analyze phase transitions in finite temperature field theory with 
multiple scalar fields. The program consists of three basic parts (see fig. [T|: 
modules for finding the tunneling solution (bubble wall profile) between dif- 
ferent vacua, a module for finding critical temperatures and phase transitions, 
and an abstract class to define specific field-theoretic models. In section |2} 
I describe the algorithms for finding bubble wall profiles for both single and 
multiple fields. Section [3] describes the algorithm for finding phase transi- 
tions, while section |4] describes how one can implement a specific model in a 
simple program. Finally, I present numerical results in section |5] and conclude 
in section ini 



To download the latest version of CosmoTransitions, visit http : / / chasm . 



ucsc . edu/ cosmotransitionSt 



2. Calculating bubble profiles 

First-order cosmological phase transitions proceed by the nucleation of 
bubbles of true vacuum out of metastable false vacuum states. The bubbles 
have both surface tension and internal pressure, so that large bubbles tend to 
expand and small bubbles tend to collapse. Critical bubbles — bubbles that 
are just large enough to avoid collapse — will drive the phase transition. 

Given a Lagrangian 

C=\{d,^){d^$)-V{$), (1) 
where is a vector of scalar fields, a critical bubble can be found by extrem- 
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specific model -^^ - » ^ generic potential .py ^ 




transitionFinder . py I ►( pathDef ormation . py 



tunnelinglD . py 



Figure 1: Overview of tlie CosmoTransitions package file structure. The 
files pathDef ormation. py and tunnelinglD. py find critical bubble profiles, 
transitionFinder .py finds the minima of finite temperature potentials as a func- 
tion of temperature and analyzes phase transitions, and generic_potential .py defines 
an abstract class that can easily be subclassed to examine a specific model. 



izing the Euclidean action 



Se = / d'^x 



(2) 



where d = 4 (3) for tunnehng at zero (finite) temperature. This quantity is 
critical for finding the nucleation rate, and thereby determining the phase 
transition temperature and whether or not the transition actually happens 
(see Refs. [311 E21 ES] for seminal work on phase transitions in quantum field 
theory). The bubble nucleation rate per unit volume is T/V = Ae""^^ at zero 
temperature and T/V = Ae~^'^^'^ at finite temperature. The prefactor A is 
quite difficult to calculate, but it has only weak temperature dependence and 
can generally be estimated on dimensional grounds (see, e.g., Ref. |31])- By 
requiring that the expectation value for one bubble to nucleate per Hubble 
volume is ~ one can show that the bubble nucleation temperature for 

weak-scale fields is given by Se/T ~ 140. 

Assuming spherical symmetry, the bubble's equations of motion are 

At finite temperature, p is simply the spatial radial coordinate and a = 2. At 
zero temperature, = r"^ — t^ and a = 3. Let 0t and (pp denote the true and 
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Figure 2: The equations of motion for a field with a potential V{4>) can be thought of as 
the equations for a particle moving in an inverted potential —V{(j)). 

false vacua, respectively. Then in order for the solution to match the field 



that bubble is non-singular at the origin. 

2.1. One- dimensional solution 

When the field has only one dimension, the solution to the critical bubble 
profile can easily be solved by the overshoot /undershoot method (see, e.g., 
Ref. |35]). Here, it is easiest to think of the problem as that of a classical 
particle moving under the influence of the inverted potential —V{(f)) plus a 
peculiar looking friction term, where (p takes on the role of a spatial coor- 
dinate and p acts as the time coordinate. The problem then is to find the 
initial placement of the particle near (px such that it rolls down the potential 
and comes to a stop at 0i? when p = oo (see figure [2]). If the particle rolls 
past (overshoots) (pp, then the initial placement was too close to 0r- If it 
doesn't have enough energy to make it to (pp (an undershoot), then the initial 
placement needs to be closer to 0t- Through trial and error, one can find 
the initial placement to arbitrary precision. 

The code presented with this paper follows the general overshoot /undershoot 
method, implemented by the class tunnelinglD.bubbleProf ile. Calcula- 
tion of thick-walled bubbles is straightforward, but thin-walled bubbles re- 
quire extra consideration. In thin-walled bubbles, the transition from cp iprp 
to ~ 0i? happens over a distance much shorter than the bubble's overall 
size. In the particle analogy, the particle sits very close to cpx for a very long 
time before quickly rolling down the potential and stopping at (pp- Extreme 
accuracy in cpQ = (p{p = 0) would be needed to reliably calculate the wall pro- 
file, since a small change in (pQ would lead to a large change in the radius of the 
bubble. Instead, I define a new variable x such that (po = (pp + e~^{(pF — (Pt), 
and use it as the initial condition to vary instead. 



at infinity, we require that 0(oo) = (pp- We also demand that ^ 



so 



p=0 
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Figure 3: A problematic potential with more than two minima. A tunneling solution from 
4>A to 4>c is only guaranteed to exist for the topmost potential (blue line). 



For small p and 



T, the equation of motion can be approximated as 



(Pel) a d(j) 
dp"^ p dp 



which has the exact solution 

0(p) -ct>T = Cp-^h [pVh) 



(4) 



(5) 



where v 



2 •■ 



and lu is the modified Bessel function of the 



first kind. For a = 2, this simplifies to (f){p) — 0t oc ^ sinh(p-\/6). These 
can be numerically inverted to find p(0), which allows one to calculate the 
approximate radius of a thin-walled bubble without performing any integra- 
tion. The integration can then start at the edge of the bubble wall, which 
both increases accuracy and decreases computation time. 

The overshoot /undershoot implementation is fairly robust, but there are 
situations in which it either does not work or is unreliable. Consider the 
potentials shown in figure |3} They each have three minima (pA < 4'b < 
with V{(j)A) > V{<Pc)- If one tries to calculate the bubble profile for tunneling 
from (pA to (pcj only the topmost potential is unproblematic. There, there 
will be a solution where the field starts near (pc, rolls over the bump at 
(pB, and ends up at (pA- In the bottommost potential, no such solution 
exists. The minimum at 0^ is the true vacuum of the theory, so the field will 
tunnel from (pA to (pB, ignoring (pc completely. The overshoot /undershoot 
method will return this solution (any initial (po > (pB will register as an 
undershoot), but it won't perform any optimization for thin- walled bubbles 
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and may have poor accuracy. In the middle case, with K(0a) > V{(Pb) > 
V{(l)c), there may or may not be a direct tunnehng solution from (pA to Re- 
using the particle analogy, it could be that any particle that starts near (pc 
and has enough momentum to get beyond (pB will necessarily overshoot (pA- 
In this case, the field can only tunnel to (pc in two steps: first by tunneling 
to (pB, and then by tunneling to (pc in a separate nucleation event. The 
algorithm in tunnelinglD.bubbleProf ile will only ever return the first of 
these transitions, and it may have poor accuracy. 

2.2. Multi- dimensional solution and path deformation 

With the solution to the one-dimensional problem in hand, we are ready to 
tackle the more challenging multi- dimensional problem. The overshoot /undershoot 
method no longer works because we required the unique topology of the one- 
dimensional case to determine whether a particular solution overshot or un- 
dershot the boundary condition at p = oo. Instead, I propose a method of 
path deformation. 

First, assume as an initial guess that the tunneling occurs on some fixed 
path in field space. That is, (p guess = (p{x), where x parametrizes the path 



dx 



1. The equations of motion split 



and for simplicity we require that 
into two parts — one parallel and one perpendicular to the path: 

dp^ p dp 
d'^d) f dx^'^ 



-2 + -„T-= in^:^)] (6) 



where V±V represents the components of the gradient of V that are perpen- 
dicular to the path. Equation [6] is the same as the one-dimensional equation 
of motion, which we can solve by the overshoot /undershoot method. If the 
path guess is correct, then the solution to equation [g] will also solve equa- 
tion [7j The trick then is to find the right path. A similar approach of path 
deformation has been proposed by Ref . , while Refs. [33 EH ES] use alter- 
nate methods to solve the equations of motion. I have previously employed 
the basic algorithm described below in the context of zero-temperature phase 
temperatures and singlet scalar dark matter models |10]. 

Using the particle analogy, we can think of the path guess as a fixed track 
on which the particle moves through the multi-dimensional space. Equation [6] 
describes the forces parallel to the track and thus determine the particle's 
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Figure 4: Path deformation in two dimensions. The normal force exerted on the starting 
path (blue straight line) pushes it in the direction of the true tunneling solution (green 
curved line). 



speed. Equation [7] does not effect the motion of the particle, but instead 
determines the normal force that the track must exert on the particle to 

keep it from falling off: = § (^g) -V±V{(f). For the right path, = 0. 
Given a starting guess, one can deform the path to the correct solution by 
continually pushing it in the direction of N (see figure |4]) . 

I implement this general method in the CosmoTransitions package using 
B-splines in the module pathDef ormation.py. Each path is written as a lin- 
ear combination of spline basis functions plus a linear component connecting 
its ends: (f){y) = f3i(f)i{y) + {(po — (f)F)y + (pp, where y parametrizes the path 
{0 <y <1, and generally g ^ 1), and 0,(0) = (^(1) = 0. This fixes the 

path's endpoints at 0o ~ 0(p = O) and (pp- Generally, only a small number 
of basis functions are needed to accurately model the path (~ 10 per field 
direction) unless it contains sharp bends or many different curves. 

Before any deformation, the algorithm first calculates the bubble profile 
along the starting path using the overshoot /undershoot method. Then it 
deforms the path in a series of steps without recalculating either the one- 
dimensional profile or ^ . At each step it calculates the normal force for 
a relatively large number of points (~ 100) along the path, rescales the 
normal force by |0r — (prl/l'^Vlmax, and moves the points in that direction 
times some small stepsize. If the one- dimensional solution is thick-walled 
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Figure 5: Example of error correction in deformation. The normal force will push wiggles 
down towards the straight line solution. But if the stepsize is too large, the wiggles will 
get reversed and amplified instead. 



(that is, (f){p = 0) is not very close 0^)? the algorithm also moves (pQ in the 
direction of averaged over the first several points (note that N{p = 0) = 
as long as the path aligns with VV^). Otherwise, 0o stays fixed at 0t- It 
then recalculates the spline coefficients by a least-squares fit to the moved 
points with the restriction that the path aligns with W at 0o when thick- 
walled. The deformation converges when the normal force is much smaller 
than forces parallel to the path. At this point, the algorithm recalculates 

^ and deforms the path again. After repeating this procedure a few times, 

';he deformation should converge after a single step and the algorithm will 
return the final tunneling solution. 

Choosing an appropriate stepsize is important. Errors in individual steps 
are generally self-correcting, but only for small stepsizes. Consider a situa- 
tion in which the correct tunneling path is a straight line, but we introduce 
errors in the deformation to add wiggles (see figure [s]). A small stepsize in de- 
formation will tend to smooth out the errors, but a large stepsize will instead 
reverse and amplify them. By checking for such reversals, the algorithm can 
keep the stepsize at the appropriate level. 



3. Exploring phase structures 

In order to determine the characteristics of a phase transition, we must 
first find where, and at what temperatures, the various phases exist. In 
theories with spontaneously broken symmetries there is at least one zero- 
temperature phase and there is generally one high-temperature symmetry- 
restoring phase. If these phases coexistence at some temperature, then there 
is likely a first-order phase transition between them. However, even in rel- 
atively simple models there can be intermediate phases (see, e.g., Ref. |41j ) 
which can lead to secondary phase transitions, or can change the quality of 
the primary transition. Therefore, it is helpful to find the location of the 
minima as a function of temperature. 



9 



By writing the potential about any point 0' as 

1 
2 



V{1 T) = a, + h{<P, - + \m,,{cP, - 00(0, - 0;.) + ■ ■ ■ , (8) 



where h = §^ and Mij = g^.g^, , one can show that there is a (nearby) 

minimum at (pmin = —M'^h. Therefore, the change in the minimum with 
respect to temperature is 

■ -M-'?l- (9) 



dT dT' 

Given a minimum at a single temperature, this allows one to find all of the 
minima of a given phase as a function of T. A singular matrix M indicates 
a rapid change in the minimum caused by either the disappearance of the 
phase or a second-order phase transition. 

Practically, it is much easier to use this equation in conjunction with a 
minimization routine than it is to integrate it by itself. The algorithm in 
transitionFinder .traceMinimum uses the Nelder-Mead downhill simplex 
method [42j to find the local minimum at particular temperatures, and then 
uses equation|9]to find how the minimum changes. The former acts as an error 
check on the latter, which allows for an adaptive stepsize in the temperature. 
Once the phase disappears, the downhill simplex method can search for a 
new phase and then trace that. In this manner, one can trace the phase 
structure of the entire theory (assuming that each phase can be found by 
minimizing the potential at one of the ends of the other phases). 

4. Structure of a simple program 

There are essentially three parts to a simple program using my code: the 
tunneling algorithms and phase tracing algorithms described above, and the 
implementation of a specific model. Much of this last task happens in the 
generic_potential class, which must be subclassed to study any particular 
theory. 

From the point of view of the finite temperature effective potential, the 
theory is completely determined by the tree-level potential and field-dependent 
mass spectrum. The generic_potential class calculates the one-loop cor- 
rections from these masses using MS renormalization 



log^ 



(10) 
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where rti and rrii are the numbers of degrees of freedom and the field-dependent 

masses of each particle species. The quantity Q is the renormahzation scale; 
c = 1/2 for gauge boson transverse modes and 3/2 for all other particles; 
and the upper (lower) sign is for bosons (fermions). The one- loop finite- 
temperature corrections are 

i 

with 

Mx) dyyHog (l T e-V^) . (12) 

The functions J^{x) are implemented in the module finiteT.py using di- 
rect integration and cubic interpolation. Of course, one can add additional 
structure (such as counter-terms) to a subclass. 

To create a fully functioning model using the Cosmo Transitions pack- 
age, one need only subclass generic.potential and overwrite three func- 
tions: the initialization function initO to specify the number of field di- 
mensions, the tree- level potential function V0(), and the mass-spectrum 
functions boson_massSq() and f ermionjnassSqO . Additionally, one should 
overwrite approxZeroTMinO to return the approximate locations of the zero 
temperature minima, especially if more than one minimum exists. Calling 
the function getPhasesO will run routines from transitionFinder .py and 
calculate the phase structure of the theory. The critical temperature (s) (the 
temperature of degenerate minima between two phases) can be found by 
calling the function calcTcTransO, while the function calcFullTrans () 
will find the amounts of supercoohng and critical bubble profiles for each 
transition. 

5. Numerical results 

5.1. Deformation 

To test the path deformation and tunneling routines, I use a simple po- 
tential given by 

V{x, y) = {x" + y^) [l.8(x - 1)^ + 0.2(y - if - 5] . (13) 

This has one local minimum at x = y = 0, and a global minimum near 
X = y = 1. For 5 <^ 1, the phases are nearly degenerate and any tunneling 
between them will be thin-walled. 



mi{(j)) 



(11) 
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X 




Figure 6: Top: the potential V{x,y) given by equation 13 for 5 = 0.02 (left) and 6 = 0.4 
(right). Black lines show successive deformations of the tunneling path, while red dashed 
lines show direct integration of equation |3] with manually chosen initial conditions. Note 
that the x and y axes are not to scale, so Vl^ does not appear perpendicular to the contour 
lines. Bottom: the bubble profiles associated with each of the above deformations, where 
4>{p) is measured along the path. 



I run the pathDef ormation . f ullTunneling class for both thin {5 = 0.02) 
and thick- walled {5 = 0.4) potentials, with results shown in figure [6] Each 
line represents 15 individual deformation steps with adaptive stepsizes on the 
order of 0.005. In this case, and in general, the thick-walled case converges 
more slowly due to the added complication of moving 0o = 0(p = O) with each 
step. The thin-walled case converges in ~ 60 steps, while the thick-walled 
case takes ~ 150 steps. 

To check the deformation solution, I numerically integrate the equations 
of motion |3j using manual trial and error to find the correct initial conditions. 
These are shown in figure |6] as the red dashed lines. The deformation algo- 
rithm gets extremely close to — but not exactly to — the integrated solution. 
In the thin- walled case, the deformed path is within 0.1% of the integrated 
solution, while in the thick-walled case the error is 0.3%. 
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Table 1: Model parameters 



mi 


m2 




V 


Q 


y\ 


y% 




120 GeV 


50 GeV 


25 GeV 


246 GeV 


246 GeV 


0.10 


0.15 


30 



5.2. Calculating transition temperatures 

To demonstrate the transition finding algorithms, I consider a theory with 
two scalar fields and a tree-level potential 

Voisu S2) = l4 i'l + (4 - v') - ^^'s,s,. (14) 

8 f ^ % 

When ix^ = 0, the theory has four identical minima located at (si,S2) = 
(iw, ±v) and (si, S2) = (±f , =Fw)- The effect of /x^ > is to lower the minima 
at {±v,±v) and raise those at (±f,=pf). The tree-level scalar mass-matrix 
is 

1 fml{3sl-v^) -/i^ 



= ^ ^"^^^::;^2 ^ ^ ^2(3,r_ j , (15) 

so that the masses at the tree-level minima are mi and m2 when /i^ = 0. I 
also add an extra bosonic degree of freedom with the field-dependent mass 

mxisuS2) = y\{sl + si) +ylsiS2. (16) 

The peculiar coupling ys is designed to lift the minima at (±f , ±v) relative 
to (±f , =pf ) at finite temperature so that there may be a phase transition 
between the two. 

I examine a model with the parameters given in table [T| where Q is the 



renormalization scale used in equation 10 and nx is the number of degrees of 



freedom assigned to the extra boson with mass mx- For this particular choice 
of parameters all tree- level minima are (meta) stable, but loop corrections 
destroy the zero-temperature minima at (±f,=Ff). I chose the parameters 
primarily to showcase multiple transitions within a single model, not for any 
physically motived purpose. The tunneling condition is Se/T = 140. 

Figure [7] shows the phase structure of the model, and figure [8] shows the 
general evolution of the potential as a function of temperature. At high 
temperatures, there is a single phase with Si = S2 = 0. When the tem- 
perature drops to 128.2 GeV, a new phase appears at (si,S2) = (43,-33) 
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GeV. The system tunnels to this phase at T = 128.1 GeV via a first-order 
phase transition (see figure [9]), and the high-temperature phase disappears 
by T = 127.6 GeV. At T = 112.2 GeV, there is a small discontinuity due 
to the non-analyticity of J_(m/T) when m = 0. At the default resolu- 
tion the algorithm registers this as a second-order phase transition, but for 
a small range of temperatures (112.14-112.22 GeV) there are actually two 
distinct phases separated by A0 ^ 2 GeV. Immediately below this the sys- 
tem can be thought of as being in the phase associated with the tree-level 
vev at (si,S2) = {+v,—v). As the system cools, the phase associated with 
(■Si, S2) = (+f , +v) appears, and by T = 75.2 GeV the two phases are degen- 
erate. Below this the former phase is metastable, but it does not transition 
until T = 54.4 GeV, which is just above the point at which it disappears. 
Thus, the final tunneling is very thick- walled (see figure |9|. 

In order to reproduce these calculations, one needs to create an instance of 
the class testModels .modell with parameters given by table [T] and call the 
class functions getPhasesO, calcTcTrans () and calcFullTrans () with 
their default parameters. Much of the plotting can be handled by the func- 
tions plot2d() and plotPhasesPhi () . 

To check that the code works with three scalar fields, I simply add an ex- 
tra field to the potential such that its minimum is always zero: Vo(si, S2, S3) = 
Vo('5i,S2) + v'^sl- I then rotate S2 and S3 by 45°. This gives an easily visu- 
alizable potential with non-trivial minima. Running the above commands 
produces the same results as in the two-dimensional model within the de- 
fault error tolerances. I successfully tested a four dimensional model in a 
similar way. In four dimensions, the transition finding routine with default 
resolution labels the highest temperature transition as second-order instead 
of weakly first-order, but this can be corrected by increasing the resolution. 
Higher dimensional models have not been tested, but the code was written 
to support an arbitrary number of scalar fields. 

6. Conclusion 

I presented the publicly available CosmoTransitions package to analyze 
cosmological phase transitions. This included algorithms to find the temperature- 
dependent phase minima, their critical temperatures, and the actual nucle- 
ation temperatures and tunneling profiles of the transitions. I introduced a 
novel method of path deformation to find the profiles, which I then demon- 
strated in simple test cases to accuracies of order ~ 0.1%. The deformation 
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Figure 7: The position of the minima as a function of temperature. 
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Figure 8: The finite-temperature effective potential V{si, S2,T) a,t T — 0, 55, 75, and 128 
GeV. Thick black lines at T = 55 and 128 GeV show the tunneling paths for the phase 
transitions. 
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Figure 9: Bubble profiles along the tunneling paths for both the high- and low-temperature 
phase transitions. The high-temperature profile is (relatively) thin-walled, while the low- 
temperature profile is very much thick- walled. 
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algorithm has been successfully tested in 2 and 3 dimensions with both thick- 
and thin-walled profiles, but it should work in any number of higher dimen- 
sions as well. 

CosmoTransitions is designed to be easily extensible, with minimal work 
needed on the part of the model builder. A new model can be created by 
subclassing the generic_potential class and specifying only it's number of 
dimensions, tree-level potential, and field-dependent particle spectrum. This 
will (hopefully) allow for the quick analysis of phase transitions in many 
extensions to the standard model. 
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